﻿using System;
using System.Collections.Generic;

namespace Allegro.Mathlib
{
    public class UniformRandom : IRandom
    {
        private int _seed;
        private Random _random;

        public UniformRandom()
            : this(DateTime.Now.GetHashCode())
        {

        }
        public UniformRandom(int seed)
        {
            _seed = seed;
            _random = new Random(seed);
        }
        #region Implementation of IRandom

        public double DrawOne()
        {
            return _random.NextDouble();
        }

        public double[] DrawMany(int n)
        {
            double[] samples = new double[n];
            for (int i = 0; i < n; i++)
            {
                samples[i] = DrawOne();
            }
            return samples;
        }

        #endregion
    }

    public class UniformRandom2 : IRandom
    {
        private const int IM1 = 2147483563;
        private const int IM2 = 2147483399;
        private const double AM = (1.0 / IM1);
        private const int IMM1 = IM1 - 1;
        private const int IA1 = 40014;
        private const int IA2 = 40692;
        private const int IQ1 = 53668;
        private const int IQ2 = 52774;
        private const int IR1 = 12211;
        private const int IR2 = 3791;
        private const int NTAB = 32;
        private const int NDIV = 1 + IMM1 / NTAB;
        private const double EPS = 1.2e-7;
        private const double RNMX = 1.0 - EPS;

        private int _idum2 = 123456789;
        private int _iy = 0;
        private int[] _iv;

        private int _seed;

        public UniformRandom2()
            : this(DateTime.Now.GetHashCode())
        {

        }
        public UniformRandom2(int seed)
        {
            if (seed == 0)
                seed = 1; // Prevent seed = 0;
            _seed = Math.Abs(seed);
            _iv = new int[NTAB];

            _Reinit();
            //_set = false;
        }

        public double DrawOne()
        {
            double temp;
            double random;
            int j, k;

            k = _seed / IQ1; // Sart here when not initializing
            _seed = IA1 * (_seed - k * IQ1) - k * IR1; // Compute seed = (IA1*seed) % IM1 without overflows by Schrage's method
            if (_seed < 0)
                _seed += IM1;
            k = _idum2 / IQ2;
            _idum2 = IA2 * (_idum2 - k * IQ2) - k * IR2; // Compute idum2 = (IA2*idum) %IM2 likewise
            if (_idum2 < 0)
                _idum2 += IM2;
            j = _iy / NDIV; // Will be in the range 0..NTAB-1
            _iy = _iv[j] - _idum2; // Here _seed is shuffled, _seed and _idum2 are combined to generate output
            _iv[j] = _seed;
            if(_iy < 1)
                _iy += IMM1;
            temp = AM*_iy;
            if(temp > RNMX)
                random = RNMX;
            else
                random = temp;
            return random;
        }

        public double[] DrawMany(int n)
        {
            double[] samples = new double[n];
            for (int i = 0; i < n; i++)
            {
                samples[i] = DrawOne();
            }
            return samples;
        }

        void _Reinit()
        {
            int j, k;

            // We don't have an extra deviate handy, so pick two uniform numers in the square extending
            // from -1 to +1 in each direction
            _idum2 = _seed;
            for (j = NTAB + 7; j >= 0; j--)
            {
                k = _seed / IQ1;
                _seed = IA1 * (_seed - k * IQ1) - k * IR1;
                if (_seed < 0)
                    _seed += IM1;
                if (j < NTAB)
                    _iv[j] = _seed;
            }
            _iy = _iv[0];
        }

    }
}
